#ifndef _WFES1_H_
#define _WFES1_H_

//#include "Equations/TwoPhaseEquation.hpp"
#include "./TwoPhaseEquation.hpp"

namespace Tuna {

  template <typename T, int Dim> class WFES1;

  template <typename T>
  struct Typeinfo <WFES1<T, 1> > {
    typedef T prec_t;
    enum { Dim = 1 };
  };

  template <typename T>
  struct Typeinfo <WFES1<T, 2> > {
    typedef T prec_t;
    enum { Dim = 2 };
  };

  template <typename T>
  struct Typeinfo <WFES1<T, 3> > {
    typedef T prec_t;
    enum { Dim = 3 };
  };
  
  template<typename Tprec, int Dim>
  class WFES1 : public TwoPhaseEquation<WFES1<Tprec, Dim> >
  {

    typedef TwoPhaseEquation<WFES1<Tprec, Dim> > TP_WFES1;

    using GeneralEquation< TP_WFES1 >::phi;
    using GeneralEquation< TP_WFES1 >::aE;
    using GeneralEquation< TP_WFES1 >::aW;
    using GeneralEquation< TP_WFES1 >::aN;
    using GeneralEquation< TP_WFES1 >::aS;
    using GeneralEquation< TP_WFES1 >::aF;
    using GeneralEquation< TP_WFES1 >::aB;
    using GeneralEquation< TP_WFES1 >::aP;
    using GeneralEquation< TP_WFES1 >::sp;
    using GeneralEquation< TP_WFES1 >::dx;
    using GeneralEquation< TP_WFES1 >::dy;
    using GeneralEquation< TP_WFES1 >::dz;
    using GeneralEquation< TP_WFES1 >::dt;
    using GeneralEquation< TP_WFES1 >::bi;
    using GeneralEquation< TP_WFES1 >::ei;
    using GeneralEquation< TP_WFES1 >::bj;
    using GeneralEquation< TP_WFES1 >::ej;
    using GeneralEquation< TP_WFES1 >::bk;
    using GeneralEquation< TP_WFES1 >::ek;
    using GeneralEquation< TP_WFES1 >::applyBoundaryConditions1D;
    using GeneralEquation< TP_WFES1 >::applyBoundaryConditions2D;
    using GeneralEquation< TP_WFES1 >::applyBoundaryConditions3D;

    using TP_WFES1::S;
    using TP_WFES1::p;
    using TP_WFES1::phi_0;
    using TP_WFES1::Srw;
    using TP_WFES1::Sro;
    using TP_WFES1::mu_w;
    using TP_WFES1::mu_o;
    using TP_WFES1::k;
    using TP_WFES1::porosity;
    using TP_WFES1::theta;
    using TP_WFES1::p_threshold;

  public:
    typedef Tprec prec_t;
    typedef typename TunaArray<prec_t, Dim >::huge ScalarField;
    
    WFES1() : TwoPhaseEquation<WFES1<prec_t, Dim > >() { }    
    ~WFES1() { };
    
    inline bool calcCoefficients1D(); 
    inline bool calcCoefficients2D() { };
    inline bool calcCoefficients3D() { };
    
    inline bool updateBoundaries1D();
    inline bool updateBoundaries2D() { };
    inline bool updateBoundaries3D() { };

    inline void printInfo() { std::cout << " WFES1 "; }
  };

/*
 *  Lineal for realtive permeability , Upwind for Sw
 */
template<typename Tprec, int Dim>
inline bool WFES1<Tprec, Dim>::calcCoefficients1D () 
{
    static prec_t Sw_e, Sw_w;
    static prec_t kdt_poro_dx_dx = k * dt / (porosity * dx * dx);
    static prec_t SrwSro1 = 1 - Srw - Sro;
    static prec_t SrwSro4 = SrwSro1 * SrwSro1 * SrwSro1 * SrwSro1;
    static prec_t SrwSro7 = SrwSro4 * SrwSro1 * SrwSro1 * SrwSro1;
    static prec_t COEF1 = kdt_poro_dx_dx / ( mu_w * SrwSro4 );
    static prec_t Sw_eSrw, Sw_eSrw4, Sw_eSrw5;
    static prec_t Sw_wSrw, Sw_wSrw4, Sw_wSrw5;
    static prec_t COEF2 = -kdt_poro_dx_dx * p_threshold / 
      (mu_w * theta * sqrt(SrwSro7));
    static prec_t aPs, aEs, aWs;

    aE = 0.0; aW = 0.0; aP = 0.0; sp = 0.0;

    for (int i =  bi; i <= ei; ++i) {      

      // Upwind
      if ( phi_0(i+1) > phi_0(i) ) Sw_e = phi_0(i+1);
      else                         Sw_e = phi_0(i);
      if ( phi_0(i-1) > phi_0(i) ) Sw_w = phi_0(i-1);
      else                         Sw_w = phi_0(i);

      Sw_eSrw  = Sw_e - Srw;
      Sw_eSrw4 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;
      Sw_eSrw5 = Sw_eSrw * Sw_eSrw4;

      Sw_wSrw  = Sw_w - Srw;
      Sw_wSrw4 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;
      Sw_wSrw5 = Sw_wSrw * Sw_wSrw4;

      aE (i) = COEF1 * Sw_eSrw4;
      aW (i) = COEF1 * Sw_wSrw4;
      aP (i) = aE (i) + aW (i);

      aEs = COEF2 * sqrt( Sw_eSrw5 );
      aWs = COEF2 * sqrt( Sw_wSrw5 );
      aPs = aEs + aWs;

      sp (i) = aPs * S(i) - aEs * S(i+1) - aWs * S(i-1);
      
    }    

/* ----- Boundary conditions ----- */


    // Neumann left side    
    Sw_wSrw  = S(bi-1) - Srw;
    Sw_wSrw5 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;
    aWs = COEF2 * sqrt(Sw_wSrw5);
    aP(bi) -= aW(bi) ;
    sp(bi) += (aWs * ( S(bi-1) - S(bi) ) + 3.53e-6 * dt / (dx * porosity)) ;
    aW(bi) = 0;

    // Update Saturation at left boundary
    /* */
    phi(bi-1) = phi(bi)
      - aW(bi) * (p(bi) - p(bi-1))// / aWs
      + 3.53e-6 * dt / (2 * porosity * dx); // * aWs)  ;
    /* */

    // Neumann right side
    Sw_eSrw  = S(ei+1) - Srw;
    Sw_eSrw5 = Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw * Sw_eSrw;
    aEs = COEF2 * sqrt(Sw_eSrw5);
    aP(ei) += aE(ei) ;
    sp(ei) += ( aEs * (S(ei+1) - S(ei)) + 2 * aE(ei+1) * p(ei+1) ) ;
    aE(ei) = 0;

    // Update Saturation at right boundary
    /* */
    phi(ei+1) = phi(ei);
    /* */
    return 0;
}


/*
 *  Update boundaries according to BC
 *
template<typename Tprec, int Dim>
inline bool WFES1<Tprec, Dim>::updateBoundaries1D () 
{
    static prec_t Sw_e, Sw_w;
    static prec_t SrwSro1 = 1 - Srw - Sro;
    static prec_t Sw_wSrw, Sw_wSrw5, 
      SrwSro7 = SrwSro1 * SrwSro1 * SrwSro1 * 
                SrwSro1 * SrwSro1 * SrwSro1 * SrwSro1;
    static prec_t COEF = - p_threshold / (mu_w * theta * sqrt(SrwSro7));
    static prec_t aPs, aEs, aWs;

    Sw_wSrw  = S(bi-1) - Srw;
    Sw_wSrw5 = Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw * Sw_wSrw;        
    aWs = COEF * sqrt( Sw_wSrw5 );

    phi_0(bi-1) = phi(bi) - aW(bi) * (p(bi) - p(bi-1)) / aWs 
                  - dx * 0.5 * 3.53e-6 / ( k * Fw);

    return 0;
}
**/


} // Tuna namespace


#endif //_WFES1_H_

















